Die DDT-Konzentrationen im Boden um eine ehemalige DDT-Fabrik in Nowshera, Pakistan können in R mit Geostatistik übersichtlich dargestellt werden. Und nicht nur das; die Konzentrationen zwischen Messpunkten werden interpoliert und eine flächendeckende Vorhersage wird in einem farbcodierten Konturenplot gezeigt.

Ablaufschema für Daten zur Auswertung mit geoR.

#— Laden des Pakets geoR —#
library(geoR)

#— Daten (t.daten) als “geodata” Objekt (d.log10conc) speichern —#
d.log10conc

#— Lagklassenbreite 40 m innerhalb des Monitoring Radius’ (t.daten$r) —#
r.bin40 pairs.min = 30, trend = ~t.daten$r)

Variogramme mit verschiedenen Lagklassenbreiten.

#— Exponentielle Schätzung mit der geeigneten Lagklassenbreite 40 m —#
r.vfit.exp #Startwert für part. Sill und Range
cov.model = “powered.exponential”,
fix.nugget = FALSE, nugget = 0, #Startwert für Nugget
fix.kappa = F, kappa = 1, #Startwert für Kappa
max.dist = 400, #max. Lagdistanz für Anpassung
weights = “equal”, #OLS
minimisation.function = “optim”, #Wahl des Optimierungslogarithmus
control = list(parscale = c(0.01,100,0.01,1)),
limits = pars.limits(kappa = c(-10,10)))

Exponentielles, sphärisches und Matern Modell passen gleich gut. Für die Vorhersage wurde das exponentielle Modell gewählt.

#— Dataframe mit Koordinaten der Gitterpunkte für Vorhersage herstellen —#
t.x t.y gitter #— Abstand der Vorhersagepunkte zur DDT-Fabrik berechnen —#
d.pred.ab

#— Kriging-vorhersage mit dem exponentiellen Modell
und lag-Klassenbreite 40 m —#
r.pred.exp krige = krige.control(obj.model = r.vfit.exp,
trend.d = ~log10(t.daten$r), trend.l = ~d.pred.ab))


#— Bestimmung des Farbkeils für Vorhersage —#
range(r.pred.exp$predict)
t.breaks t.col.scale

#— Kriging Vorhersage —#
image(r.pred.exp, locations = gitter, values = r.pred.exp$predict,
breaks = t.breaks, x.leg = c(-180,180), y.leg = c(-500,-450), col = t.col.scale, xlab = “x-Coordinate [m]”,
ylab = “y-Coordinate [m]”)


#— vertikale & horizontale Linie, gemessene Konzentrationen, Kontourlinien der Konzentrationen, Beschriftung —#
abline(v = 0, lty = 2)
abline(h = 0, lty = 2)
points(t.daten$x, t.daten$y, cex = 0.2*t.daten$conc)
r.cont.z contour(t.x, t.y, z = r.cont.z, add = T, labcex = 0.8)
title( main = “DDT Kriging Prediction”)
text(-280, -425, “log DDT [mg/kg]”, cex = 0.8)